home *** CD-ROM | disk | FTP | other *** search
/ Info-Mac 4 / Info_Mac IV CD-ROM (Pacific HiTech Inc.)(August 1994).iso / Science / RLaB / rlib / lu.r < prev    next >
Text File  |  1994-04-25  |  2KB  |  92 lines
  1. //-------------------------------------------------------------------//
  2. //  Syntax:    lu ( A )
  3.  
  4. //  Description:
  5.  
  6. //  Lu computes the "LU" factors of the input MATRIX. The input MATRIX
  7. //  must be square. lu() returns l, u, and pvt in a LIST. The
  8. //  factorization has the form: 
  9.  
  10. //  A = p*l*u            where A is the input matrix
  11.  
  12. //  The function lu is a user-function that utilizes the built-in
  13. //  function factor. Factor does the factorization using the LAPACK
  14. //  functions DGETRF/ZGETRF, and DGECON/ZGECON.
  15.  
  16. //  Example:
  17.  
  18. //  > a = [1,2,3;4,5,6;7,8,0]
  19. //   a =
  20. //   matrix columns 1 thru 3
  21. //             1           2           3
  22. //             4           5           6
  23. //             7           8           0
  24. //  > f = lu(a)
  25. //    l            pvt            u            
  26. //  > f.p
  27. //   matrix columns 1 thru 3
  28. //             0           0           1
  29. //             1           0           0
  30. //             0           1           0
  31. //  > f.l
  32. //   matrix columns 1 thru 3
  33. //             1           0           0
  34. //         0.143           1           0
  35. //         0.571         0.5           1
  36. //  > f.u
  37. //   matrix columns 1 thru 3
  38. //             7           8           0
  39. //             0       0.857           3
  40. //             0           0         4.5
  41. //  > f.p*f.l*f.u
  42. //   matrix columns 1 thru 3
  43. //             1           2           3
  44. //             4           5           6
  45. //             7           8           0
  46.  
  47. //  See Also: backsub, factor, inv, solve
  48. //-------------------------------------------------------------------//
  49.  
  50. static (swap);
  51.  
  52. lu = function ( A )
  53. {
  54.   local (i, l, u, pvt, x)
  55.  
  56.   if (A.nr != A.nc) { error ("lu() requires square A"); }
  57.  
  58.   x = factor (A);    // Do the factorization
  59.  
  60.   //
  61.   // Now create l, u, and pvt from lu and pvt.
  62.   //
  63.  
  64.   l = tril (x.lu, -1) + eye (size (x.lu));
  65.   u = triu (x.lu);
  66.   pvt = eye (size (x.lu));
  67.  
  68.   //
  69.   // Now re-arange the columns of pvt
  70.   //
  71.  
  72.   for (i in 1:max (size (x.lu)))
  73.   {
  74.     pvt = pvt[ ; swap (1:pvt.nc, i, x.pvt[i]) ];
  75.   }
  76.   return << l = l; u = u; pvt = pvt >>;
  77. };
  78.  
  79. //
  80. //  In vector V, swap elements I, J
  81. //
  82.  
  83. swap = function ( V, I, J )
  84. {
  85.   local (v, tmp);
  86.   v = V;
  87.   tmp = v[I];
  88.   v[I] = v[J];
  89.   v[J] = tmp;
  90.   return v;
  91. };
  92.